home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / lisp / gcl-1.000 / gcl-1 / gcl-1.0 / lsp / numlib.lsp < prev    next >
Encoding:
Lisp/Scheme  |  1994-05-07  |  6.5 KB  |  219 lines
  1. ;; Copyright (C) 1994 M. Hagiya, W. Schelter, T. Yuasa
  2.  
  3. ;; This file is part of GNU Common Lisp, herein referred to as GCL
  4. ;;
  5. ;; GCL is free software; you can redistribute it and/or modify it under
  6. ;;  the terms of the GNU LIBRARY GENERAL PUBLIC LICENSE as published by
  7. ;; the Free Software Foundation; either version 2, or (at your option)
  8. ;; any later version.
  9. ;; 
  10. ;; GCL is distributed in the hope that it will be useful, but WITHOUT
  11. ;; ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. ;; FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public 
  13. ;; License for more details.
  14. ;; 
  15. ;; You should have received a copy of the GNU Library General Public License 
  16. ;; along with GCL; see the file COPYING.  If not, write to the Free Software
  17. ;; Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19.  
  20. ;;;;    numlib.lsp
  21. ;;;;
  22. ;;;;                           number routines
  23.  
  24.  
  25. (in-package 'lisp)
  26. (export
  27.  '(isqrt abs phase signum cis asin acos sinh cosh tanh
  28.    asinh acosh atanh
  29.    rational rationalize
  30.    ffloor fround ftruncate fceiling
  31.    lognand lognor logandc1 logandc2 logorc1 logorc2
  32.    lognot logtest
  33.    byte byte-size byte-position
  34.    ldb ldb-test mask-field dpb deposit-field
  35.    ))
  36.  
  37.  
  38. (in-package 'system)
  39.  
  40.  
  41. (proclaim '(optimize (safety 2) (space 3)))
  42.  
  43.  
  44. (defconstant imag-one #C(0.0d0 1.0d0))
  45.  
  46.  
  47. (defun isqrt (i)
  48.        (unless (and (integerp i) (>= i 0))
  49.                (error "~S is not a non-negative integer." i))
  50.        (if (zerop i)
  51.            0
  52.            (let ((n (integer-length i)))
  53.                 (do ((x (ash 1 (ceiling n 2)))
  54.                      (y))
  55.                     (nil)
  56.                     (setq y (floor i x))
  57.                     (when (<= x y)
  58.                           (return x))
  59.                     (setq x (floor (+ x y) 2))))))
  60.  
  61. (defun abs (x)
  62.        (if (complexp x)
  63.            (sqrt (+ (* (realpart x) (realpart x))
  64.                     (* (imagpart x) (imagpart x))))
  65.            (if (minusp x)
  66.                (- x)
  67.                x)))
  68.  
  69. (defun phase (x)
  70.        (atan (imagpart x) (realpart x)))
  71.  
  72. (defun signum (x) (if (zerop x) x (/ x (abs x))))
  73.  
  74. (defun cis (x) (exp (* imag-one x)))
  75.  
  76. (defun asin (x)
  77.        (let ((c (- (* imag-one
  78.                       (log (+ (* imag-one x)
  79.                               (sqrt (- 1.0d0 (* x x)))))))))
  80.             (if (or (and (not (complexp x))
  81.               (<= x 1.0d0)
  82.                (>= x -1.0d0)
  83.                 )
  84.                 (zerop (imagpart c)))
  85.                 (realpart c)
  86.                 c)))
  87.  
  88. (defun acos (x)
  89.        (let ((c (- (* imag-one
  90.                       (log (+ x (* imag-one
  91.                                    (sqrt (- 1.0d0 (* x x))))))))))
  92.             (if (or (and (not (complexp x))
  93.               (<= x 1.0d0)
  94.                (>= x -1.0d0)
  95.                 )
  96.                 (zerop (imagpart c)))
  97.                 (realpart c)
  98.                 c)))
  99.  
  100.  
  101. (defun sinh (x) (/ (- (exp x) (exp (- x))) 2.0d0))
  102. (defun cosh (x) (/ (+ (exp x) (exp (- x))) 2.0d0))
  103. (defun tanh (x) (/ (sinh x) (cosh x)))
  104.  
  105. (defun asinh (x) (log (+ x (sqrt (+ 1.0d0 (* x x))))))
  106. (defun acosh (x)
  107.        (log (+ x
  108.                (* (1+ x)
  109.                     (sqrt (/ (1- x) (1+ x)))))))
  110. (defun atanh (x)
  111.        (when (or (= x 1.0d0) (= x -1.0d0))
  112.              (error "The argument, ~s, is a logarithmic singularity.~
  113.                     ~%Don't be foolish, GLS."
  114.                     x))
  115.        (log (/ (1+ x) (sqrt (- 1.0d0 (* x x))))))
  116.  
  117.  
  118. (defun rational (x)
  119.   (etypecase x
  120.     (float      
  121.       (multiple-value-bind (i e s) (integer-decode-float x)
  122.                (if (>= s 0)
  123.                    (* i (expt (float-radix x) e))
  124.                  (- (* i (expt (float-radix x) e))))))
  125.     (rational x)))
  126.  
  127.  
  128. (setf (symbol-function 'rationalize) (symbol-function 'rational))
  129.  
  130. ;; although the following is correct code in that it approximates the
  131. ;; x to within eps, it does not preserve (eql (float (rationalize x) x) x)
  132. ;; since the test for eql is more strict than the float-epsilon
  133.  
  134. ;;; Rationalize originally by Skef Wholey.
  135. ;;; Obtained from Daniel L. Weinreb.
  136. ;(defun rationalize (x)
  137. ;  (typecase x
  138. ;    (rational x)
  139. ;    (short-float (rationalize-float x short-float-epsilon 1.0s0))
  140. ;    (long-float (rationalize-float x long-float-epsilon 1.0d0))
  141. ;    (otherwise (error "~S is neither rational nor float." x))))
  142. ;
  143. ;(defun rationalize-float (x eps one)
  144. ;  (cond ((minusp x) (- (rationalize (- x))))
  145. ;        ((zerop x) 0)
  146. ;        (t (let ((y ())
  147. ;                 (a ()))
  148. ;             (do ((xx x (setq y (/ one
  149. ;                                   (- xx (float a x)))))
  150. ;                  (num (setq a (truncate x))
  151. ;                       (+ (* (setq a (truncate y)) num) onum))
  152. ;                  (den 1 (+ (* a den) oden))
  153. ;                  (onum 1 num)
  154. ;                  (oden 0 den))
  155. ;                 ((and (not (zerop den))
  156. ;                       (not (> (abs (/ (- x (/ (float num x)
  157. ;                                               (float den x)))
  158. ;                                       x))
  159. ;                               eps)))
  160. ;                  (/ num den)))))))
  161.  
  162.  
  163. (defun ffloor (x &optional (y 1.0s0))
  164.        (multiple-value-bind (i r) (floor (float x) (float y))
  165.         (values (float i r) r)))
  166.  
  167. (defun fceiling (x &optional (y 1.0s0))
  168.        (multiple-value-bind (i r) (ceiling (float x) (float y))
  169.         (values (float i r) r)))
  170.  
  171. (defun ftruncate (x &optional (y 1.0s0))
  172.        (multiple-value-bind (i r) (truncate (float x) (float y))
  173.         (values (float i r) r)))
  174.  
  175. (defun fround (x &optional (y 1.0s0))
  176.        (multiple-value-bind (i r) (round (float x) (float y))
  177.         (values (float i r) r)))
  178.  
  179.  
  180. (defun lognand (x y) (boole boole-nand x y))
  181. (defun lognor (x y) (boole boole-nor x y))
  182. (defun logandc1 (x y) (boole boole-andc1 x y))
  183. (defun logandc2 (x y) (boole boole-andc2 x y))
  184. (defun logorc1 (x y) (boole boole-orc1 x y))
  185. (defun logorc2 (x y) (boole boole-orc2 x y))
  186.  
  187. (defun lognot (x) (logxor -1 x))
  188. (defun logtest (x y) (not (zerop (logand x y))))
  189.  
  190.  
  191. (defun byte (size position)
  192.   (cons size position))
  193.  
  194. (defun byte-size (bytespec)
  195.   (car bytespec))
  196.  
  197. (defun byte-position (bytespec)
  198.   (cdr bytespec))
  199.  
  200. (defun ldb (bytespec integer)
  201.   (logandc2 (ash integer (- (byte-position bytespec)))
  202.             (- (ash 1 (byte-size bytespec)))))
  203.  
  204. (defun ldb-test (bytespec integer)
  205.   (not (zerop (ldb bytespec integer))))
  206.  
  207. (defun mask-field (bytespec integer)
  208.   (ash (ldb bytespec integer) (byte-position bytespec)))
  209.  
  210. (defun dpb (newbyte bytespec integer)
  211.   (logxor integer
  212.           (mask-field bytespec integer)
  213.           (ash (logandc2 newbyte
  214.                          (- (ash 1 (byte-size bytespec))))
  215.                (byte-position bytespec))))
  216.  
  217. (defun deposit-field (newbyte bytespec integer)
  218.   (dpb (ash newbyte (- (byte-position bytespec))) bytespec integer))
  219.